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Abstract: A common approach for generating an 
anisottopic mesh is the M-uniform mesh approach where 
an adaptive mesh is generated as a uniform one in the 
metric specified by a given tensor M. A key component 
is the determination of an appropriate metric which is of- 
ten based on some type of Hessian recovery. Recently, 
the use of a global hierarchical basis error estimator was 
proposed for the development of an anisotropic met- 
ric tensor for the adaptive finite element solution. This 
study discusses the use of this method for a selection of 
different applications. Numerical results show that the 
method performs well and is comparable with existing 
metric tensors based on Hessian recovery. Also, it can 
provide even better adaptation to the solution if applied 
to problems with gradient jumps and steep boundary 
layers. For the Poisson problem in a domain with a cor- 
ner singularity, the new method provides meshes that are 
fully comparable to the theoretically optimal meshes. 

Keywords: mesh adaptation, anisotropic mesh, finite 
element, a posteriori estimate, hierarchical basis, vari- 
ational problem, anisotropic diffusion. 

MSC 2010: 65N30, 65N50. 
1 Introduction 

A common approach for generating an anisotropic mesh 
is the M-uniform mesh approach based on generation 
of a quasi-uniform mesh in the metric space defined 
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by a symmetric and strictly positive definite metric ten- 
sor M. A scalar metric tensor will lead to an isotropic 
mesh while a full metric tensor will generally result in 
an anisotropic mesh. In this sense, the mesh generation 
procedure is the same for both isotropic and anisotropic 
mesh generation. A key component of the approach is 
the determination of an appropriate metric often based 
on some type of error estimates. 

Typically, the appropriate metric tensor depends on 
the Hessian of the exact solution of the underlying prob- 
lem, which is often unavailable in practical computation, 
thus requiring the recovery of an approximate Hessian 
from the computed solution. A number of recovery tech- 
niques are used for this purpose, for example the gradi- 
ent recovery technique by Zienkiewicz and Zhu Il35ll36l . 
the technique based on the variational formulation by 
Dolejsi 1 18 1, or the quadratic least squares fitting (QLS) 
proposed by Zhang and Naga ll34l . Generally speak- 
ing, Hessian recovery methods work well when exact 
nodal function values are provided (e.g. interpolation 
problems), but unfortunately they do not provide an ac- 
curate recovery when applied to linear finite element ap- 
proximations on non-uniform meshes, as pointed out by 
the author in [29]. Recently, conditions for asymptot- 
ically exact gradient and convergent Hessian recovery 
from a hierarchical basis error estimator have been given 
by Ovall IIBTI . His result is based on superconvergence 
results by Bank and Xu f8][9l, which require the mesh 
to be uniform or almost uniform: assumptions which are 
usually violated by adaptive meshes. 

Hence, a convergence of adaptive algorithms based 
explicitly on the Hessian recovery cannot be proved in a 



1 



direct way, even if tiieir application is successful in prac- 
tical computations [ 1 8l |26l |30l |33l . This explains the 
recent interest in anisotropic adaptation strategies based 
on some type of a posteriori error estimates. For ex- 
ample, Cao et al. fT3l studied two a posteriori error 
estimation strategies for computing scalar monitor func- 
tions for use in adaptive mesh movement; Apel et al. 0] 
investigated a number of a posteriori strategies for com- 
puting error gradients used for directional refinement; 
and Agouzal et al. HI El El and Agouzal and Vassilevski 
H proposed a new method for computing metric ten- 
sors to minimize the interpolation error provided that an 
edge-based error estimate is given. 

Recently, Huang et al. [24 1 presented a mesh adapta- 
tion method based on hierarchical basis error estimates 
(HBEE). The new framework is developed for the linear 
finite element solution of a boundary value problem of a 
second-order elliptic partial differential equation (PDE), 
but it is quite general and can easily be adopted to other 
problems. A key idea in the new approach is the use of 
the globally defined HBEE for the reliable directional 
information: globally defined error estimators have the 
advantage that they contain more directional informa- 
tion of the solution; error estimation based on solving lo- 
cal error problems, despite its success in isotropic mesh 
adaptation, do not contain enough directional informa- 
tion, which is global in nature; moreover, Dobrowolski 
et al. ifTTl have pointed out that local error estimates can 
be inaccurate on anisotropic meshes. A brief description 
of the method is provided in Sect. |2] 

The objective of this article is to study the application 
of this new anisotropic adaptation approach to different 
problems. 

The first example deals with a boundary value prob- 
lem of a second-order elliptic PDE. Numerical results in 
Il24l have already shown that the new mesh adaptation 
approach can be a successful alternative to Hessian re- 
covery in mesh adaptation for a boundary value problem 
of a second-order elliptic PDE. In this paper, we would 
Uke to investigate another interesting question, namely 
whether the proposed algorithm is able to provide op- 
timal mesh adaptation for a problem where the optimal 
mesh design is known theoretically. For this purpose, a 
Poisson problem in a domain with a corner singularity is 
studied, which requires proper mesh adaptation in order 
to obtain an accurate finite element solution. Theoret- 
ically optimal meshes are known for this problem and 
are based on some a priori information of the domain 



and the solution. Section [3] investigates if the proposed 
method is able to generate adaptive meshes that are com- 
parable to the optimal ones without any a priori infor- 
mation on the exact solution. A brief overview over the 
considered mathematical problem, the properties of the 
optimal meshes and the numerical results for the new 
method are provided. 

Section |4] presents an anisotropic metric tensor for 
general variational problems developed by Huang et al. 
f25\ using the HBEE and the underlying variational 
formulation and gives a numerical example for a non- 
quadratic variational problem. The metric tensor is com- 
pletely a posteriori: it is based solely on the residual, 
edge jumps, and the a posteriori error estimate. 

The third example is an anisotropic diffusion prob- 
lem. The exact solution of this problem satisfies the 
maximum principle and it is desirable for the numerical 
solution to fulfill its discrete counterpart: the discrete 
maximum principle (DMP). Recently, Li and Huang 
li30ll developed an anisotropic metric tensor based on the 
anisotropic non-obtuse angle condition, which provides 
both mesh adaptation and DMP satisfaction for the nu- 
merical solution: the mesh alignment is determined by 
the main diffusion drag direction, i.e. by the underly- 
ing PDE, and the Hessian of the exact solution deter- 
mines the optimal mesh density. In Sect. |5] the Hessian 
of the exact solution is replaced with the Hessian of the 
hierarchical error estimator to obtain a new, completely 
a posteriori, metric tensor accounting for both DMP sat- 
isfaction and mesh adaptation. 

Concluding remarks on the numerical examples and 
some key components of the hierarchical basis error es- 
timator are given in Sect. |6] 

2 Anisotropic mesh adaptation based on 
hierarchical basis error estimator 

Consider the solution of a variational problem: find u G 
V such that 

a{u,v)=f{v) VvGV (P) 

where V is an appropriate Hilbert space of functions 
over a domain H G M^, a(-, •) is a bilinear form defined 
onV xV, and F(-) is a continuous linear functional on 
V. 

For a given simplicial mesh the linear finite ele- 
ment approximation m/, of u is the solution of the cor- 
responding variational problem in a finite dimensional 
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subspace Vh CV of piecewise linear functions over 
find Mft G Vh such that 

a{uh,Vh) = f{vh) Vv/, e Vh- (Ph) 

For the adaptive finite element solution, (and thus V/,) 
is generated according to a given quantity of interest, for 
example the error of the solution in a chosen norm. 

This study follows the M-uniform mesh approach 
ll23l which generates an adaptive mesh as a uniform 
mesh in the metric specified by a symmetric and strictly 
positive definite tensor M = M{x). Such a mesh is called 
an M-uniform mesh. Once a metric tensor M has been 
chosen, a sequence of mesh and corresponding finite el- 
ement approximation are generated in an iterative fash- 
ion. 

Let {^1^'^} (/ = 0, 1,...) be an affine family of sim- 
plicial meshes on Q. with the corresponding space vj^''' 
of continuous, piecewise hnear functions over An 
adaptive algorithm starts with an initial mesh ^3^^^°' . On 
every mesh ^^j-'^ the variational problem {Pf,) with vj^''' 

is solved and the obtained approximation u^^ is used to 
compute a new adaptive mesh for the next iteration step. 
The new mesh =^3/,''^^' is generated as a M-uniform mesh 

with a metric tensor Mj^ defined in terms of m['^ . This 
yields the sequence 

The conformity of the mesh to the input metric is char- 
acterized by the mesh quality measure Qmesh introduced 
in (241. Qmesh > 1 for any given mesh but Q,„esh = 1 if 
and only if the underlying mesh is M-uniform (see ll24l 
Sect. 4.1] for more details). In the presented numeri- 
cal examples, the mesh adaptation process is repeated 
until the mesh is M-uniform within a given tolerance: 
Qmesh < 1 + £ with a chosen tolerance £ (£ = 0.1 is used 
in the numerical experiments throughout the paper). 

Mesh generation software BAMG is used to generate 
new adaptive meshes {bidimensional anisotropic mesh 
generator developed by F. Hecht Il22l ). 

2.1 Adaptation based on a posteriori error 
estimates 

Typically, the metric tensor M/, depends on the Hessian 
of the exact solution of the underlying problem [20, .23]. 
As mentioned in the introduction, it is not possible to 



obtain an accurate Hessian recovery from a linear finite 
element solution in general |29|, so there is no way to 
prove a convergence of an adaptive algorithm based on 
the Hessian recovery in a direct way, even if its applica- 
tion is successful in practical computations |[T8ll26ll30l . 

An alternative approach developed in fl^ employs an 
a posteriori error estimator for defining and computing 
M/,. The brief idea is as follows. 

Assume that an error estimate Zh is reliable in the 
sense that 

< C\\zh\\ ■ (1) 
for a given norm || • || and that it has the property 

n^zh = (2) 

for some interpolation operator 11^. Then the finite ele- 
ment approximation error is bounded by the (explicitly 
computable) interpolation error of the error estimate Zh, 
viz., 

<C'||z/,|| =C||z/, -nftz/,||. (3) 

Now, it is known from the interpolation theory fZTl that 
the interpolation error for a given function v can be 
bounded by a term depending on the triangulation 
and derivatives of v, i.e., 

||v-n/,v|| <c^(^ft,v), 

where C is a constant independent of and v. There- 
fore, we can rewrite Q as 

< CS'{3^h,Zh)- 

In other words, up to a constant, the solution error is 
bounded by the interpolation error of the error estimate. 
Thus, the metric tensor M/, can be constructed to min- 
imize the interpolation error of the zn and does not de- 
pend on the Hessian of the exact solution. 

2.2 Hierarchical basis a posteriori error estimate 

One possibility to achieve the property (|2]l is to use the 
hierarchical basis error estimator. The general frame- 
work can be found among others in the work of Bank 
and Smith |7| or Deuflhard et al. 1 16,1 . The approach is 
briefly explained as follows. 
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Let eh = u — M;, be the error of the Unear finite element 
solution M/, G Vft. Then for all v G V we have 

a{eh,v)= f{v)-a{uh,v). (E) 

Let Vi, = Vh(B Wi, be a space of piecewise quadratic 
functions, where Wh is the linear span of the quadratic 
edge bubble functions (a quadratic edge bubble function 
is defined as a product of the two linear nodal basis func- 
tions corresponding to the edge endpoints). The error 
estimate Zh is defined as the solution of the approximate 
error problem: find Zh S Wh such that 

a{zh,Wh)=f{wh)-a{uf„Wh) Vw/, GW),. 

The estimate Zh can be viewed as a projection of the true 
error onto the subspace W/, and relies on two assump- 
tions: the strengthened CBS-inequality 

\a{vh,Wh)\ < Y\\\vh\\\ \\\wh\\\ with 7< 1 

for all V/, G V/, and w/, G Wh and the saturation assumption 

III" — w^lll < jSIIIm — with j3 < 1, 

i.e. the assumption that the quadratic finite element ap- 
proximation Uq is more accurate than the lineai" approx- 
imation M/,. 

Now, if Ylf, is defined as the vertex-based, piecewise 
Unear Lagrange interpolation then zn satisfies the condi- 
tion Q since the edge bubble functions vanish at ver- 
tices. Then, if assumption ([T]) holds, the finite element 
approximation error can be controlled by minimizing the 
interpolation error of z/,, i.e., the right-hand side in Q. 

Note that this definition of the error estimate is global 
and the cost of its exact solution is comparable to the 
cost of the quadratic finite element approximation. To 
avoid the expensive exact solution in numerical compu- 
tation, only a few sweeps of the symmetric Gauss-Seidel 
iteration are employed for the resulting linear system 
(until the relative difference of the old and the new ap- 
proximations is under a given relative tolerance). In nu- 
merical experiments, three to ten sweeps were enough to 
achieve the relative tolerance of 1% and it proved to be 
fully sufficient for the purpose of mesh adaptation. The 
computational cost of such approximation is compara- 
ble to the cost of the Hessian recovery: in the tests, the 
computation of HBEE was in average about two times 
slower than the Hessian recovery method. Since both 
methods require only a fraction of the overall computa- 
tional cost, using the fast HBEE approximation for mesh 
adaptation instead of Hessian recovery will not signifi- 
cantly increase the overall computational time. 



3 Poisson problem in a domain with a 
corner singularity 

Consider the Dirichlet problem for the Poisson equation 

Uu = f in a, 
]u = g ondD., 

with a = {x : r < I, < d < Tl/X} for some 
^<A<l,/ = and g = r'^ sin in usual polar co- 
ordinates (r, 0). The solution of Q or, more precise, the 
solution of the corresponding variational problem 
given by 

u{r,e) = r^ sin(A0). 

The gradient of the solution u has a singular be- 
haviour near the corner (0,0) of the domain Q. and 
u G H^+^^^{Q.) with arbitrarily small e L28il . 

3.1 Optimal mesh grading 

For a uniform or quasi-uniform mesh, the poor regular- 
ity of u leads to the suboptimal rate of convergence of 
the linear finite element method: 

\\it-Uh\\H<{Q.) ^Ch'^^^ ^ (5) 
ll"-"/illL2(n) <C/^^^"''- (6) 

A number of techniques were developed to regain the 
optimal convergence order by using specially adapted fi- 
nite element spaces, see 15, Sect. 4.2] and the references 
therein for more details. In particular, it is possible to 
regain the optimal convergence rate by using the stan- 
dard linear finite element space with proper mesh grad- 
ing around the corner of the domain. 

The basic idea is described among others in ll28l and 
is as follows. The linear finite element error is bounded 
by the linear interpolation error 

\u — Uh\Hi(Q,)^\l^ — '^hU\Hi{Q.)^C f^K\'^\H^{K) ■ 

Ke.% 

Thus, in order to obtain the optimal convergence rate, 
the mesh size hk has to be balanced with the size of 
\^\h-(k) ■ Now, if mesh elements in a neighbourhood 
around the corner (0,0) have the size 

hK = Ch/f^, 
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where h is the mesh size far away from the comer and 
vk is the distance of an element K to the corner, then the 
hnear finite element error is 

^ Ch, 

h-Uh\\L2(^ci) <Ch^: 

provided that < A lH |28l. Note, that such meshes 
still have A/^ ~ h^^ number of elements, thus significantly 
increasing the accuracy of the approximation for a given 
number of mesh elements. 

In this example, we follow the approach described in 
the Sect. |2] and employ the metric tensor developed in 
IJ?! for a boundary value problem of a second-order el- 
liptic PDE. In two dimensions and for the L^-norm of 
the error, the metric tensor Mub based on the KB error 
estimator is given element-wise as 
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MHB,K = A&t\l+ — \HK{Zh. 



I+ — \HK{Zh 

«/1 



where Hnizh) is the Hessian of the quadratic hierarchi- 
cal basis error estimate z.h on element K and «/, is a reg- 
ularization parameter to ensure that the metric tensor is 
strictly positive definite. Also, «/; can be used as a adap- 
tation intensity control: if «/, — )• °°, the mesh becomes 
uniform; if «/, — )• 0, the mesh becomes more adaptive. 
Usually, ah is chosen so that about half of the mesh ele- 
ments are concentrated in regions where det(M) is large 
(see ||24|| for more details on the choice of a/,). 

Note, that Mhb is completely a posteriori and does 
not require any a priori information on the solution. 

3.2 Numerical example 

Consider the Dirichlet boundary problem Q with A = 
4/7. In two dimensions, the number of elements N of a. 
quasi-uniform mesh is proportional to h^^, so there is a 
factor of —1/2 in the convergence rate if it is given in 
terms of the number of mesh elements. Thus, accord- 
ing to ([5]) and Q, the expected orders of convergence in 
terms of the number of mesh elements for quasi-uniform 
meshes should be at most —2/7 (i.e. —A /2) for the 
and —4/7 (i.e. —A) for the norms of the error. For the 
optimally graded meshes we should expect orders —0.5 
and — 1 , respectively. 

Examples of uniform and adaptive meshes as well as 
close-up views near the corner are given in Figs.[T]and[2] 
For the adaptive mesh, the concentration of mesh points 




(a) Quasi-uniform, 1214 triangles. (b) M//b, 1225 triangles. 

Figure 1: Corner singularity and optimal mesh grading: 
mesh examples. 
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(a) Quasi-uniform: 165936trian- (b) M//b> 166 125 triangles, 
gles. 

Figure 2: Corner singularity and optimal mesh grading: 
close-ups at (0,0). 



near the singularity can be clearly observed (Figs, lb 



and 2b I. We can also observe the gradual change of 



mesh elements size in the close-up of the adaptive mesh 



(Fig. 2b I ; in particular, along the boundary edges from 
the domain corner to the outer boundary. The mesh 
grading is rather moderate, but it turns out to be enough 
to achieve the optimal convergence order, as can be ob- 
served in the corresponding convergence plot (Fig. [3]). 

Figure [5] shows the and norms of the global 
linear finite element solution error on uniform and adap- 
tive meshes against the number of mesh elements. As 
expected, the convergence order of the finite element so- 
lution on uniform meshes is only —2/7 and —4/7 for the 
and the L? norms of the error, respectively. Contrary, 
the convergence order of the error of the finite element 
solution obtained by means of the metric tensor Mhb is 
-1/2 and -1. 

Thus, the considered method is able to achieve the 
optimal mesh grading in this example relying on only on 
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h\ uniform — 
L^, uniform — ^— 
h\Mhb -h- 




10^ lO"* 10^ 

Number of elements N = 0(ln'^) 



Figure 3: Corner singularity and optimal mesh grading: 
global error for quasi-uniform and adaptive 
meshes against the number of elements. 



the a posteriori information provided by the hierarchical 
basis error estimator. 



4 Variational problems 

While it has attracted considerable attention from many 
researchers and been successfully applied to the numer- 
ical solution of PDEs, anisotropic mesh adaptation has 
rarely been employed for variational problems, espe- 
cially when combined with a posteriori error estimates. 
Recently, Huang and Li [,26] developed a metric ten- 
sor for the adaptive finite element solution of variational 
problems. In the anisotropic case, it is semi-a posteri- 
ori: it involves residual and edge jumps, both dependent 
on the computed solution, and the Hessian of the exact 
solution. In [25], this result was improved to provide 
a metric tensor for variational problems based on the 
HBEE and the underlying variational formulation. The 
new metric tensor is a posteriori in the sense that it is 
based solely on residual, edge jumps, and HBEE. This 
is in contrast to most previous work where M depends on 
the Hessian of the exact solution and is semi-a posteriori 
or completely a priori; e.g., see llTTl[T2l[T4ll2Tll23ll26ll . 

4.1 General variational problem and the 
anisotropic metric tensor 

Consider a general functional of the form 

I[v] = f F{x,v,yv)dx, Vv G Vg 



where F {■,■,■) is a given smooth function, Q. C U.'^ 
(d = 1,2,3) is the physical domain and Vg is a properly 
selected set of functions satisfying the Dirichlet bound- 
ary condition 

v{x)=g{x) yxGdO. 

for a given function g. 

The corresponding variational problem is to find a 
minimizer u £Vg such that 

I[u]=mml[v]. 

veVg 

A necessary condition for m to be a minimizer is that the 
first variation of the functional vanishes. This leads to 
the Galerkin formulation 

5l[u,v] = / Vm) v + Fvj,(x,m,Vm) • Vv)(3fx = 

Ja 

(7) 

for all V G Vb, where Vb = Vg with g = and F„ and Fy^ 
are the partial derivatives of F with respect to u and Vu, 
respectively. 

Given a triangulation ^ for Q. and the associated lin- 
ear finite element space Vg^ C Vg, the finite element 
solution M/, can be found by solving the corresponding 
Galerkin formulation: find m/, G Vgji such that 

/ {Fu {x, Uh ,Vuh)vh+ Fsju {x,Uh, Vuh ) • Vv/, )dx = 
Jo. 

for all Vfi G Vb,/i. 

The a posteriori metric tensor Mhb,k for general vari- 
ational problems developed in |[25]| for the error mea- 
sured in the //^-semi-norm is given element-wise by 



M, 



HB,K 



1 + 



1 



\K 



1/2 



1 1/2 



yedK 

xA&t[I+ — \HK{Zh 




I+ — \HKiZh) 
CCf, 



with the residual 

rH{x)=F,{x)-V-Fyu{x) xeK e 
and the edge jump 

Rh{x) = (Fv„(x)-?iy)|^ + (Fv„(x)-?iy)|A'/ 



6 




Figure 4: Variational problem: numerical solution. 

on 7 G d£^h\dD. (/?/, = if 7 G dD.). As in the previous 
section, ai, is a regularization parameter to ensure that 
the metric tensor is strictly positive definite (see [24] for 
more details on the choice of «/,). 

Note, that dl[u,v] in ^ is linear in v but is nonlin- 
ear in u in general. Thus, a modification of the error 
problem ^Ef^ for Zh in Sect. |2.2| is required. For this pur- 
pose, denote by ah{ui,;-, •) a bilinear form resulting from 
a linearization of •] about m/, with respect to the first 
argument. The error estimate Zh is then defined as the 
solution of the approximate linear error problem: find 
Zh £ such that 

ah{uh,Zh,Wh) = -dl[uh,Wh] 

for allwA e W/, Q. 

4.2 Numerical example 

Consider an anisotropic variational problem defined by 
the non-quadratic functional 



I[u]=J (^l + |VMp)^^V 1000m? 



dx 



with = (0, 1) X (0, 1) and the boundary condition 

{u = I onx = Oorx=l, 
u = 2 onj = Oorj = l. 



This example is discussed in 11251 126II and is originally 
taken from ifTOl ; the analytical solution is not available, 
but a computed solution in Fig. |4] shows that the mesh 
adaptation challenge for this example is the resolution 
of the sharp boundary layers near x = and x= I. 




(a) Metric based on residual, edge jumps, and Hessian recovery: 
1160 triangles, maximum aspect ratio 15. 




(b) Metric based on residual, edge jumps, and HBEE: 1143 trian- 
gles, maximum aspect ratio 51. 

Figure 5: Variational problem: adaptive meshes and 
close-up views at (0,0). 



Adaptive meshes obtained by means of Hessian re- 
cover){^ and HBEE are given in Fig. [s] The both 
methods have correct mesh concentration and provide 
good alignment with the boundary layers. Anisotropic 
meshes are comparable, although mesh elements near 
the boundary layer in the HBEE-based adaptive mesh 
have a larger aspect ratic|^ than elements of the mesh 
obtained by means of the Hessian recovery. This could 
be due to the smoothing nature of the Hessian recovery: 
usually, it operates on a larger patch, thus introducing an 
additional smoothing effect, which affects the grading of 
the elements' size and orientation. The global hierarchi- 
cal basis error estimator does not have this handicap and, 
in this example, the mesh obtained by means of HBEE 
has a better refinement in the orthogonal direction along 
the steep boundary layers. 



' Quadratic least squares fitting Hessian recovery 1341 . which seems 
to be the most robust and reliable Hessian recovery method 1291 

^In this paper, aspect ratio is defined as the longest edge divided by 
the shortest altitude. An equilateral triangle has an aspect ratio 



of2/V3^ 1.15. 
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5 Anisotropic diffusion and the DMP 

Anisotropic diffusion problems arise in various areas of 
science and engineering, for example image process- 
ing, plasma physics, or petroleum engineering. Stan- 
dard numerical methods can produce spurious oscilla- 
tions when they are used to solve these problems. A 
common approach to avoid this difficulty is to design a 
proper numerical scheme or a mesh so that the numeri- 
cal solution satisfies the discrete counterpart (DMP) of 
the maximum principle satisfied by the continuous so- 
lution. A well known condition for the DMP satisfac- 
tion by the linear finite element solution of isotropic 
diffusion problems is the non-obtuse angle condition 
that requires the dihedral angles of mesh elements to be 
non-obtuse |[T5l . In |[30l . a generalization of the con- 
dition, the so-called anisotropic non-obtuse angle con- 
dition, was introduced for the finite element solution of 
heterogeneous anisotropic diffusion problems. The new 
condition is essentially the same as the existing one ex- 
cept that the dihedral angles are measured in a metric de- 
pending on the diffusion matrix of the underlying prob- 
lem, i.e. the mesh is aligned with the major diffusion 
directions. Based on the new condition, a metric tensor 
for anisotropic mesh adaptation was developed, which 
combines the satisfaction of the DMP with mesh adap- 
tivity: in two dimensions and for the error measured in 
the //'-semi-norm, it is given element- wise by 



Mdmp,k 



1 



1 + —Bh,k] det(Dif)2D^\ (8) 



where 



B 



H.K 



detl 



1 

1^ 



\H(u)\\rdx 



(9) 

and ah is a regularization parameter to ensure that the 
metric tensor is strictly positive definite (see [24l for 
more details on the choice of a/,). The diffusion ma- 
trix D in ([8]) provides the correct mesh alignment and 
thus the satisfaction of the DMP whereas the Hessian 
of the exact solution in Q provides adaptivity with the 
solution profile. 

As noted in Sect. [2} if the finite element solution error 
can be bounded by the interpolation error of the hierar- 
chical basis error estimator, the Hessian of the exact so- 
lution in ([9]) can be replaced by the Hessian of the HBEE 





M = * " " X 

(a) Boundary conditions. (b) Numerical solution. 

Figure 6: Anisotropic diffusion example. 



Zh- 



BHB,K = det{OK) 



K Jk 



}K\H{zh)\fdx. 

(10) 

With this choice, the mesh will still satisfy the DMP be- 
cause the mesh is still aligned with the major diffusion 
directions, but this time the mesh density is determined 
by the a posteriori error estimate z/,. 

5.1 Numerical example 

Consider the EVP discussed in |[30l : 

fV-(DVM)=/ inQ., 
\u = g on do., 



with 

/ = 0, a=[0,l]2\ 



4 5 

9'9 



n 2 



g = Oonr„„,, g = 2on T,-, 



where Tout and r,„ are the outer and inner boundaries of 



D., respectively (Fig. 6a). The diffusion matrix is given 
by 



cos G — sin 
sin 6 cos 6 



1000 
1 



cos 6 sind 
— sin 6 cos 6 



with the angle of the primary diffusion direction 

6 = 7rsinA;cos3' 

(the primary diffusion direction is parallel to the first 
eigenvector of D). 

This example satisfies the maximum principle and the 
solution stays between and 2 and has sharp jumps near 
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Figure 7: Anisotropic diffusion: meshes and contour Figure 8: Anisotropic diffusion: meshes and con- 
tour plots of the numerical solution for the 



plots of the numerical solution for the (a) 



isotropic and (b) HBEE-based (Mhb) metric 
tensors. 



DMP-compliant (a) Hessian recovery-based 
(Mdmp+h) and [(b)] HBEE-based (Mdmp+hb) 
metric tensors. 



the inner boundary. The analytical solution is not avail- 



able, but the numerical solution is provided in Fig. 6b 



The goal is to produce a numerical solution which also 
satisfies DMP, i.e. stays between and 2, and has a good 
adaptation. 

To emphasize the compliance with the DMP, the met- 
ric tensors Mjjmp+h and Mdmp+hb based on (|9]l and 
( fTO] ), respectively, are also compared to a uniform mesh 
and the anisotropic metric tensor Mhb based solely on 
the HBEE (Sect. [3]). Figures [7] and [8] show meshes and 
solution contours. 

No overshoots in the finite element solutions are ob- 
served for all cases, but Fig. [7] shows that undershoots 
and unphysical minima occur in the solutions obtained 
with the uniform mesh (minw/, ^ —0.059) and Mhb 
(minw/, w —0.0032). As expected, no undershoots can 
be observed for Momp+h and Mdmp+hb (Fig-[8])- 

As for MoMP and Momp+hb, the solution contours for 
both metric tensors are almost undistinguishable and, al- 
though not as smooth, still quite comparable to the one 



obtained with Mhb (cf. Figs. [8] and 7b I, thus providing 
a good adaptation to the sharp solution jump near the 
interior boundary (cf. the somewhat smeared solution 
jump for the isotropic mesh in FigjTaj right). The mesh 
computed by means of HBEE is fully comparable with 
the mesh obtained using Hessian recovery. However, the 
maximum aspect ratio of the mesh obtained by means of 
the HBEE is, again, slightly larger. 

6 Concluding remarks 

Numerical results show that a global HBEE can be a 
successful alternative to Hessian recovery in mesh adap- 
tation: the new method performs successfully and is 
quite comparable with the commonly used Hessian re- 
covery method. A fast approximate solution is as fast 
as Hessian recovery and proved to be sufficient to pro- 
vide enough directional information for the purpose of 
the mesh adaptation. Moreover, as already observed in 
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Il24l . it could be of advantage for problems with discon- 
tinuities, because Hessian recovery could result in un- 
necessarily high mesh density for such problems. The 
global HBEE seems to be less affected by this issue 
and, depending on the underlying problem, can pro- 
vide a more appropriate mesh adaptation. Also, Hes- 
sian recovery can also cause a light mesh smoothing: 
meshes obtained by means of Hessian recovery-based 
method in Sect. |4] have a smaller maximum aspect ratio 
than meshes obtained with HBEE and therefore seem to 
be slightly worse in terms of adaptation with the steep 
boundary layers. For the comer singularity problem, the 
discussed mesh adaptation method proved to be able to 
capture the required information and to provide the op- 
timal mesh grading without any a priori information on 
the solution. 

One of the key components of the method is the reli- 
ability of the error estimator on anisotropic meshes: er- 
ror estimation with hierarchical bases is usually based 
on the saturation assumption, which basically states that 
quadratic approximations provide finer information on 
the solution than linear ones. Existing results on its va- 
Udity require bounds on the elements' aspect ratio |19| . 
It is still unclear if similar results can be achieved for 
general adaptive meshes, but numerical results suggest 
that aspect ratio bounds are not necessary if the mesh 
is properly aligned. Moreover, it seems that good mesh 
adaptation does not require an accurate Hessian recov- 
ery or an accurate error estimator, but rather some addi- 
tional information of global nature, although it is still 
unclear which information exactly is necessary. This 
question is definitely of interest for future investigations. 

Acknowledgements. The author is very grateful to the 
anonymous referees for their valuable comments and 
suggestions. 
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